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Abstract 

We present GW calculations of molecules, ordered and disordered solids and interfaces, 
which employ an efficient contour deformation technique for frequency integration, and do 
not require the explicit evaluation of virtual electronic states, nor the inversion of dielectric 
matrices. We also present a parallel implementation of the algorithm which takes advantage of 
separable expressions of both the single particle Green's function and the screened Coulomb 
interaction. The method can be used starting from density functional theory calculations per¬ 
formed with semi-local or hybrid functionals. We applied tire newly developed technique to 
GW calculations of systems of unprecedented size, including water/semiconductor interfaces 
with thousands of electrons. 


1 Introduction 


The accurate description of the excited state properties of electrons plays an important role in many 
fields of chemistry, physics, and materials scienceFor example, the interpretation and prediction 
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of photoemission and optoelectronic spectra of molecules and solids rely on the ability to com¬ 
pute transitions between occupied and virtual electronic slates from first principles, as well as their 
lifetimes 2 . 

In particular, in the growing field of materials for energy conversion processes - including solar 
energy conversion by the photovoltaic effect and solar to fuel generation by water photocatalysis 
- it has become key to develop predictive tools to investigate the excited state properties of nanos¬ 
tructures and nanocomposites and of complex interfaces ' - '. The latter include solid/solid and 
solid/liquid interfaces, e.g. between a semiconductor or insulator and water or an electrolyte 6-10 . 
In the last three decades. Density Functional Theory (DFT) has been widely used to compute 
structural and electronic properties of solids and molecules 11-15 . Although successful in describ¬ 
ing ground state and thermodynamic properties, and in many ab initio molecular dynamics stud¬ 
ies 11 ’ 17 . DFT with both semi-local and hybrid functionals has failed to accurately describe excited 
state properties of several materials ls . However, hybrid functionals have brought great improve¬ 
ment for properties computed with semi-local ones. e.g. for defects in semiconductor and ox¬ 
ides 1 ' 1 22 . In particular hybrid functionals with admixing parameters computed self-consistently 
have shown good performance m reproducing experimental band gaps and dielectric constants of 
broad classes of systems 23 . In the case of the electronic properties of surfaces, interfaces (and 
hence nanostructures), the use of hybrid functionals has in many instances not been satisfactory. 
Indeed calculations with hybrid functionals yield results for electronic levels that often depend on 
the mixing parameter chosen for the Hartree-Fock exchange: such parameter is system dependent 
and there is no known functional yielding satisfactory results for the electronic properties of inter¬ 
faces composed of materials with substantially different dielectric properties, as different as those 
of. e.g. water (£*> = 1.78) 21 and Si <£„ = 11.9) 2 ' or water and transition metal oxides of interest 
for photoelectrodes (£, = 5-7) 2h . 

The use of many body perturbation (MBPT) starting from DFT single particle states has proven 
accurate for several classes of systems 27-36 and it appears to be a promising framework to describe 
complex nanocomposites and heterogeneous interfaces. MBPT for the calculations of photoemis- 
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sion spectra in the GW approximation 17 , and of optical spectra by solving approximate forms of 
the Bethe Salpeter Equation I BSE)- 11 ' is in principle of general applicability: however its use has 
been greatly limited by computational difficulties in solving the Dyson's equation and the BSE for 
realistic systems. 

Recently we proposed a method to compute quasi particle energies within the GqWq approxima¬ 
tion <i.e. the non-selfconsistent GW approximation) that does not require the explicit calculation of 
virtual electronic states, nor the inversion of large dielectric matrices 39,40 . In addition the method 
does not use plasmon pole models but instead frequency integrations are explicitly performed and 
there is one single parameter that controls the accuracy of the computed energies, i.e. the number 
of eigenvectors and eigenvalues used in the spectral decomposition of the dielectric matrix at zero 
frequency. The method was successfully used to compute the electronic properties of water 41 and 
aqueous solutions 42 and of heterogeneous solids 5 , including crystalline and amorphous samples 40 . 
However the original method contained some numerical approximations in the calculations of the 
head and wings of the polarizability matrix: most importantly the correlation self-energy was com¬ 
puted on the imaginary axis and obtained in real space by analytic continuation. Finally, although 
exhibiting excellent scalability, the method was not yet applied to systems with thousands of elec¬ 
trons. e.g. to realistic interfaces, due to the lack of parallelization in its original implementation. 

In this paper we solved all of the problems listed above, by (i) eliminating numerical approxima¬ 
tions in the calculation of the polarizability; (ii) avoiding the use of an analytic continuation and 
using efficient contour deformation techniques: (iii) providing a parallel implementation of the 
algorithm based on separable forms of both the single particle Green function and the screened 
Coulomb interaction. The method presented here can be used starting from DFT orbitals and ener¬ 
gies obtained both with semi-local and hybrid functionals. We applied our technique to the calcula¬ 
tion of the electronic properties of systems of unprecedented size, including water/semiconductor 
interfaces with more than one thousand electrons. These calculations allow one to accurately study 
states at heterogeneous interfaces and to define an electronic thickness of solid/liquid interfaces 
using MBPT. 
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The rest of ihe paper is organized as follows. Sec. Section 2 describes the GqWq methodology that 
we implemented in a computational package called West. Sec. Section 3 presents results for the 
ionization potentials of closed and open shell molecules and for the electronic structure of crys¬ 
talline. amorphous and liquid systems, aimed at verifying and validating the algorithm and code 
West against previous calculations and measurements. Sec. Section 4 presents the study of the 
electronic properties of finite and extended large systems, i.e. nanocrystals and solid/liquid inter¬ 
faces. of interest to photovoltaic and photocatalysis applications, respectively. Our conclusions are 
reported in Sec. Section 5. 
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2 Method 


Wilhin DFT. the /i-th single particle state Ynka and energy of a system of interacting electrons 
is obtained by solving the Kohn-Sham (KS) equation 11-1 

fifes I ty'nko) = e «k<7 I Ynko) < 1) 

where — T + P/ on + +- is the KS Hamiltonian. T is the kinetic energy operator, and 

% m , and are the ionic. Hartree. and exchange-correlation potential operators, respectively. 
The indexes k and a label a wavevector within the first Brillouin /.one (BZ) and spin polarization, 
respectively. Here we consider collinear spin. i.e. decoupled up and down spins. 

In a fashion similar to Eq. (Eq. ( 1 )) one may obtain quasiparticle (QP) states y^k<r ai, d QP energies 
by solving the equation: 



where the QP Hamiltonian is formally obtained by replacing, in Eq. (Eq. (1)), the exchange- 
correlation potential operator with the electron self-energy operator = jCMVT: G° is the in¬ 
teracting one-particle Green's function. IV is the screened Coulomb interaction and T is the vertex 
operator 28 43 . All quantities entering the definition of the self-energy are interdependent and can 
be obtained self-consistently adopting the scheme suggested by L. Hedin 44-46 . In the GW ap¬ 
proximation. T is set equal to the identity, which yields the following expression for the electron 
self-energy 47 : 

+oo ^ 

E°(r,r';G>) = i j ‘^-G a (r.r f :(o + (o')W R p A {r.v , -.a/) , (3) 

where Wrp\ is the screened Coulomb interaction obtained in the random phase approximation 
(RPA). Due to the non-locality and the frequency dependence of the electron self-energy, a self- 
consistent solution of Eq. (Eq. (2» is computationally very demanding also for relatively small 
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systems, containing tens of electrons, and usually one evaluates QP energies perturbatively: 




(4) 


£,*, + (i«,ic»|£‘ r (£,^)|yi,k 0 )-(v'.tol'^lv,*,i. (5) 


We note that E^ a enters both the left and right hand side of Eq. (Eq. (3)). whose solution is 
usually obtained recursively, e.g. with root-finding algorithms such as the secant method. The use 
of Eq. (Eq. (3)) to evaluate QP energies from KS states and of the corresponding KS wavefunctions 
is known as the GoVPo approximation. 

Within GqW o. using the Lehmann's representation, the Green’s function is: 


G£ s (r.r';*>) = -£/ 

"HZ 


<Jk V.*<J(v)v; lka (r , ) 

(2 Jty' e„ka -to- /»fsign(e nko - e F ) 


( 6 ) 


where t] is a small positive quantity and Ep is the Fermi energy. In Eq. (Eq. (6)) we used the 
subscript KS to indicate that the Green's function is evaluated using the KS orbitals obtained by 
solving Eq. (Eq. (1)). 

In Ref. (39.40| an algorithm was introduced to compute the self-energy matrix elements of Eq. (Eq. (3)) 
without the need to evaluate explicitly empty (virtual) electronic states, by using a technique called 
projective eigendecomposition of the dielectric screening (PDEP). A diagram of the method is re¬ 
ported in Fig. Figure 1. After KS single particle orbitals and energies are obtained using semilocal 
or hybrid functionals, the screened Coulomb interaction is computed using a basis set built from 
the eigenpotentials of the static dielectric matrix at zero frequency. In this way Wpp A entering 
Eq. (Eq. (3)) is expressed in a separable form, similar to that of G% s in Eq. (Eq. (6)>. In the follow¬ 
ing sections we describe in detail all the steps outlined in Fig. Figure 1. The separable form of W/fp A 
is given in Sec. Section 2.1. Calculation of the polarizability and the spectral decomposition of the 
dielectric matrix are described in Sec. Section 2.2 and Section 2.3. respectively. Matrix elements of 
Gks and Wrpa are then obtained without the explicit use of empty electronic states and simultane¬ 
ously at several frequencies by using a deflated Lanczos technique, described in Sec. Section 2.4. 
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Finally the frequency integration is carried out by introducing a contour deformation method, as 
described in Sec. Section 2.5. The use of tire analytic continuation used in the original method of 
Ref. 139.401 is thus avoided. 
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Figure I: (Color online) Schematic representation of the steps involved in the calculations of 
quasiparticle (QP) energies, within the GoWfa approximation, using the method proposed in this 
work. The KS energies (e,) and occupied orbitals <t^.) computed at the DFT level are input to the 
PDEP algorithm, which is used to iteratively diagonalize the static dielectric matrix (£ 1 ) at zero 
frequency. The set of eigenvectors {<?»} constitutes the basis set used to compute both G and VP at 
timte frequencies with die Lanczos algorithm. The frequency integration of Eq. <Eq. (3)) is carried 
out using the contour deformation technique. The frequency dependent matrix elements of the 
electron self-energy are thus obtained and introduced in Eq. (Eq. (3)) to compute the QP energies 

Ef. 


2.1 Separable form of the screened Coulomb interaction 

In order to solve equation Eq. (Eq. (3)) and obtain QP energies, one needs to compute the ma¬ 
trix elements of die electron self-energy between KS states, which in the GqVVq approximation is 
given by Eq. (Eq. (3)). The Green’s function may be expressed in a fully separable form using its 
Lehmann’s representation. Eq. |Eq. (6». However is not trivially separable: it is given as the 
sum of two terms: 

b’m(r.r':w) = v(r.r') + VP /l (r.r':u) (7) 
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where v(r.r') = -p-pf is the hare Coulomb interaction ls and W p is a nonlocal and frequency de¬ 
pendent function. Using Eq. (Eq. (7)). we write the self-energy as the sum of two contributions 
T. a — + ££• where only the latter depends on the frequency: 


Ej(r.r-) 


‘ I ^C^(r.r’:6> + w')v(r:r'> 

— oe 

,v< “ r dk 

-£ / ~2~ TV'«ka(r)v(r.r'>vr,; k « T (r') 
n=x bz “ 


( 8 ) 

(9) 


Kc is the number of occupied states with spin a: Ey is usually called exchange self-energy 
because it is formally equivalent to the Fock exact exchange operator* 9 : 


Ec(rV;ffl) = i / -r—G a K . 


( 10 ) 


E£ is referred to as correlation self-energy. Using Eq.s (Eq. (7))-(Eq. (10)) the QP Hamiltonian of 
Eq. (Eq. (2)) may be expressed as: 


V(“) = f + *um + *&, + £?(») (11) 

where *Sf is the Hartree-Fock potential operator. The ionic potential V,,*, is treated within the 
pseudopotential approach'". 

In this work we express VV' /( in a separable form by adopting the projective dielectric eigendecom- 
position (PDEP) technique, proposed in Ref. [51-52]. and we use a plane wave basis set to express 
the single particle wave functions and charge density, within periodic boundary conditions: 

VWr) = <’ ,k VWD - Ec Bk<T (G)e /(k+G) r (12) 

<; 


where G is a reciprocal lattice vector. fniur(G) = ^ u„^ a {r)e l< ’ r and Q is the unit cell 
volume. In Eq. (Eq. (12)) all reciprocal lattice vectors such that \ |k + G 1 < E culw f c are included 



in ihe summation. Using a plane wave description also for IV’,, we have 


»Wr. r; (0) = / t^tt I <* ,,q+G) - r [ V GG' + ^ G ,(q: to)] e~ l ^ r ' < 13) 

HZ { ’ 


where =- <| lx ( f .; (^ is the Kronecker delta) and 




Xvc.lq :a>) 

«, + G|q + G' 


(14) 


In Eq. (Eq. (14)) we have introduced the symmetrized reducible polarizability X- related to the 
symmetrized inverse dielectric matrix £ 1 by the relation: 




(15) 


The symmetrized form X of the polarizability X * s 

v/4; ur /A no 1 

XGG ' ~ jq + Gf^lq + G 7 ! ‘ 


( 16 ) 


The reducible polarizability X is related to the irreducible polarizability Xo by the Dyson’s equa¬ 
tion. which within the RPA reads: 


XgG , -XgG'+ Y *GG, v GiG 2 *G 2 G' 
G|,Gj 


(17) 


or in terms of symmetrized polarizabilites: 




(18) 


Within a plane wave representation each quantity in Eq. (Eq. (18)) is a matrix of dimension i\~ H , 
and in principle x c an be obtained from x" v > a simple linear algebra operations. In practice, 
the matrices of Eq. (Eq. (18)) may contain millions of rows and columns for realistic systems: 
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for example for a silicon nanocrystal with 465 atoms, placed in a cubic box of edge 90bohr. 1.5 
million plane waves are needed in the expansion of Eq. (Eq. (12)) with E cu iwfc ~ 25 Ry. It is thus 
important to find alternative representations of X and reduce the number of elements to compute. 
One could think of a straightforward spectral decomposition: 

N e**p 

*GG'(q;a>)= £ ^(q + G;o»)^(q;fi})^;(q + G;6>) (19) 

1=1 

where 0 , and A, are the eigenvectors and eignvalues of £. respectively. Unfortunately this strategy 
is still too demanding from a computational standpoint, as it implies finding eigenvectors and 
eigenvalues at multiple frequencies. 

A computationally more tractable representation may be obtained using the spectral decomposition 
of Xo at (0 — 0. As apparent from Eq. (Eq. (18)), eigenvectors of X are also eigenvectors of £°; 
the latter is easier to iteratively diagonalize than X • and the frequency dependency may be dealt 
with iterative techniques, stalling from the solution at (0 = 0. as discussed in Sec. Section 2.4. 
Therefore we proceed by solving the secular equation for x" only at (0 — 0. generating what we 
call the PDEP basis set { ify) : i — 1. p }. which is used throughout this work to represent the 
polarizability X- 

XgcA*<») = £ ^(q + GIA^q^fi^q + G); (20) 

1 = 1 . 

J=i 

here A ,j (q: to) is a matrix of dimension N^j ep . In general N p j ep <g Af ;m .. 51 '" leading to substantial 
computational savings. 1 " The N^e,, functions 0, may be computed by solving the Stemheimer 
equation 54 , without explicitly evaluating empty (virtual) electronic states. In addition. turns 
out to be the only parameter that controls the accuracy of the expansion in Eq. (Eq. (20)). The 
details of the derivation of the PDEP basis set are given in Sec. Section 2.3. We note that alternative 
basis sets, based on the concepts related to maximally localized Wannier functions, have been 
proposed in the literature to reduce the dimensionality of the polarizability matrices 55 . 
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By defining t),<q + G) = 


1+<i| ’ 


we formally obtain the desired separable form for W p : 



( 21 ) 


The scaling operation used to define t>, is divergent in the long wavelength limit <q —> 0) and for 
G — 0. However such divergence can be integrated yielding: 

W r (r.r':io) = Z((0) + - £ 0 , <r)A(r') , (22) 

“ i=i. 

>=* 

where 

^q-tl 

In Eq. (Eq. (23)) the integration is evaluated on the region ftq-o of the BZ enclosing the T-point 
(i.e. q — 0). 56 

In the q —> 0 limit, we can now write the matrix elements of Ic using: i) the separable form of 
Wp of Eq. (Eq. (22)) and ii) the expression of Gas- given in Eq. (Eq. (6)). in terms of projector 
operators: 


n/. 


where 

(25) 

t*t a = l5i Vnka) <V6,ko-| and = Ll= N o + 1 IV'-.ko) (YWl are the projector operator over the 
occupied and unoccupied manyfold of states belonging to k-point k and spin a. respectively 57 . 
Hence we have: 


(V' (1 k«r|5:r(«)IV'i,k 0 ) =^,ka(0» + S«k«T(w)+Ck O («) +fW<6>), (26) 
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where A nka and <?, lko (B llka and D, lko ) are contributions to the correlation self-energy originating 
from occupied (empty ) states: 


-*-oo ^ 

= i I ^-2(6)') (ftfarl /?*<)& (fl) + *>' + irj) /*° I(28) 
4~°° t V 

Gte(») = 0/ £ £ (29) 

—oo ^ • 

/=> 

+«° ^ 

A:k«r(» = 71 r^r L A > / ( 6 >')<&k«r| (*>+©' + in)«&*) (30) 

;=• 

We have defined 0,' ko (r) - y/„ko (r) $} (r>. The quantities A„ k(7 . B„ ka , C llk<7 and D nko entering 
Eq. (Eq. (26)) are now in a form where iterative techniques (see Sec. Section 2.4) can be applied to 
obtain the matrix elements of the correlation self-energy. Moreover, because of the completeness 
of energy eigenstates (f^ a = 1 — P^ c ), we may compute all quantities in Eq.s (Eq. (27))-(Eq. (30)) 
considering only occupied states. The integration over the frequency domain will be discussed in 
Sec. (Section 2.5). 

2.2 Polarizability within the random phase approximation 

Here we discuss how to compute the polarizability x from x° within the RPA. in the long wave¬ 
length limit (q — > 0). without explicitly evaluating electronic empty states. The Fourier components 
of the symmetrized irreducible polarizability x" are given by the Adler-Wiser expression' 8 * 9 . 
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which contains an explicit summation over unoccupied states: 


JlW -q + G q+G 7 - 


a n=\m=NZ,~\ h7 


1 


+ 


Iq + GHq + G' 

I 


L £ mka ~ ^«ik~ qo + W — / 1] £ m ka _ fuk-qo — 0) ~ IIJ ] 


(31) 


where the matrix element 


Pinnko(qiG) — (Vn.-kcl 4 *’ ' | V'nk-qo ) 
is often referred to as oscillator strength: it has the following properties: 


(32) 


P„,„k<T«q-G = 0)| q _ 0 = (33) 

^qp.m l k O (q.G = 0)| q _ 0 = /(V'mkd r|v6, k(T ). (34) 


Following Ref. (60.61], we partition the polarizability of Eq. (Section 2.2) into head (G - G' - 0). 
wings <G - 0.G' ^ 0 or G f- 0.G' - 0) and body <G f- 0 and G' f 0) elements. The q -» 0 
limit of the body, which we call is analytic, while the limits of the head and wings are 

non-analytic. i.e. they depend on the Cartesian direction along which the limit is performed. The 
long wavelength limits of the head, body and wings of the polarizability matrix are summarized in 
Table Table 1. Using the PDEP basis set we obtain: 


Table 1: The long wavelength limit (q —» 0) of the head, wing and body elements of the polar¬ 
izability Xixi'i®) are S* ven > n ihe second and third columns: Uc,p(lo) — and 

F a p(o) — 4 Ke'Xoo( a ) are evaluated using Eq. (Eq. (32)) and yield the linear and quadratic 
terms in the Taylor expansion of around q = 0. respectively. 


£r.G'(q "+ 0:fl> ) 

g' = o 

G'^O 

G = 0 

Lap l /aF a piO)qp/lf 

-iLa<laUacAto)/<l 

G #0 


Bgg 1 


13 








(35) 


u aJ ((o) = £<W(6>)$,(G') 

G' 


&',(<») = I ftWBGcAcotfjiG') 

GC 


We can now express all ihe quantities in Table Table 1 without any explicit summation over empty 
(virtual) states: 


F c,{i (W) = 4.TC- Y, £ / pJjT : [^5 ( f -'ko - 0> + llj) + + 0) + iTj)] ^ 1£,fk< 


(6>) = 4ff < r££ / (&U - <» + if?) + Ofo(e„w<! + «+'*?)] /*° 

^ w HZ 


t \ a 

( a) - 4,TC 2 £ £ J ^3 (Kto i /*° (f..k«7 - (0 + itj ) + 0£ s (f„k<T + © + *■*?)] 'Ml 


( 2*) 3 


/C. , ,/ k 

$>(<») = 4,t<’ 2 £ £ / I /?® [6% s (f„ k< r - w + irj) + O a KS {e nka + (o +it])] f* c g' k 


Note that the greek letters a and fi identify Cartesian directions, while the roman letters i and 
; label the eigevectors of x" at 0 ) - 0. i.e. the elements of the PDEP basis set. We have also 
defined the auxiliary functions ^(r) - V'i»ka(r)6(r) and <=,? k(T (r) - F* a r a % ka ). Within pe¬ 
riodic boundary conditions the position operator is ill-defined and ^£ a (r) is obtained by solving 


the linear system 


(»*s - £*.) 14k.) - ft" Iw*.) 


where the commutator of the KS Hamiltonian with the position operator includes the contribu¬ 
tion of the non-local part of the pseudopotential 60,6 Once x" is obtained, x is computed using 
Eq.(Eq.(18)). 
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The quantities required to evaluate Eq. (Eq. (26)) are the following 


h2; 

1-Alto) [ tl q 4 tt e 1 
/ ( 2 *)' |q|- 

**q—<1 


( 42 ) 


*») = |1 - BloJI-'Bla) + ^^[1 -B(<9)l->/i(<a)[l - B(ffl)]- 1 (43) 

with k(io) - 1 - /(to)-Tr{ji(to)(l - B(to)] -1 }: f{(0) = 7 L 0 / r aa<to) and the matrix elements 
of ni,(io) — Y.u Uia(®)Uaj((0)- In Eq- (Eq. (43)) the bold symbols denote matrices of dimension 
A ffx/rp- In order to compute the matrix elements of the correlation self-energy. Eq. (Eq. (26)). we 
need to evaluate E(to) and A (to), namely the head and body of the X operator. These are easily 
obtained via linear algebra operations from F a p(6)). U a j(o), U, a (o)) and B,j((0). 

By replacing explicit summations over unoccupied states with projection operations. Eq.s (Eq. (37))- 
(Eq. (40)) may be evaluated using linear equation solvers and (owing to the completeness of the en¬ 
ergy eigenstates) the calculation of polarizabilities is carried out without the explicit evaluation of 
the virtual states. In a similar fashion one obtains the auxiliary functions £L(0 ^ Eq. (Eq. (41)) 
and the PDEP basis set as described in Sec. Section 2.3. We note that other approaches were 
developed in the literature 61 66 to improve the efficiency of Go Wo calculations by avoiding the 
calculation of virtual states, or by limiting the number of virtual states to be computed. However 
such approaches did not make use of the spectral decomposition of the irreducible polarizability 
to obtain the reducible polarizability, but instead inverted explicity large matrices. Specifically in 
Reining el al. ( " the Sternheimer equation was used to obtain the irreducible polarizability without 
virtual states and then a plasmon pole model was adopted to compute the dielectric response as 
a function of frequency. In Giustino el al M the Sternheimer equation was used as well to obtain 
the irreducible polarizability without computing virtual states: the polarizability matrix was then 
inverted numerically and either a plasmon pole model or a a Pade expansion were used to treat 
the frequency dependence. In our approach we avoided large matrix inversion by using the PDEP 
basis set to express all polarizability matrices. Finally, we note that an additional advantage of our 
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approach is that Eq.s (Eq. (37))-(Eq. (40)) may be computed using a deflated Lanczos algorithm 
for multiple values of the frequency, as discussed in Sec. Section 2.4. A Lanczos algorithm was 
also used by Soininen et al 69 to iteratively include local field effects in RPA Hamiltonians and 
avoid explicit inversion of large matrices. However the authors of Ref. (69] computed explicitly 
virtual states. 

2.3 Projective dielectric eigenpotential (PDEP) basis set 

We now describe in detail how to obtain the PDEP basis set {|0;) : i — 1 .N r j ell ) introduced in 
Eq. (Eq. (20)); each function 0, is computed by the iterative diagonalization procedure, summa¬ 
rized in Fig. Figure 2. the procedure is initiated by building an orthonormal set of N p j ep basis 
vectors, e.g. with random components. Then Stemheimer equations are solved in parallel, 
where the perturbation is given by the /-th basis set vector 0,(r). In particular, given a perturbation 
the linear variation Ayr' |k(T ) of the occupied eigenstates of the unperturbed system 
may be evaluated using the Stemheimer equation 54 : 


Ws - w) |A«U = -ffvr ttto) . (441 


Eq. (Eq. (44)) may be iteratively solved using e.g. preconditioned conjugate-gradient methods. 
The linear variation of the density due to the /-th perturbation is obtained within density functional 
perturbation theory 70,71 (DFPT) as 


A/M.ri 


% *17 

= I £ / AVlaiDYnkoir) +C.C.] . 

17 "= X HZ ' } 


(45) 


The matrix elements of the irreducible polarizability in the space spanned by 0, are given by: 


X" = 4;re 2 (0, An/) 


(46) 
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where is computed using Eq.s (Eq. <44))-<Eq. (45)) and assuming that l' |KI, (G) — <p,(G). 

The matrix then diagonalized to obtain new Npj ep basis vectors 0 ,. and the procedure is it¬ 
erated using e.g. a Davidson algorithm (See Fig Figure 3). We note that at each iteration, all 
Stemheimer problems are independent from each other, thus offering the opportunity to carry out 
embarrassingly parallel calculations. A description of the parallel operations and data layout will 
be given elsewhere As a result, the algorithm shows a good scalability up to 524288 cores (see 
Fig. Figure 4). 

As an example we show in Fig. Figure 5 the eigenvalues of the matrix obtained with the PDEP 
algorithm for the water, silane, benzene and sodium chloride molecules, using KS Hamiltoni¬ 
ans with different exchange-correlation functionals. The choice of the functional only affects the 
most screened eigenpotentials. whereas the eigenvalues A, corresponding to the least screened ones 
rapidly approach 51 zero with a decay similar to that predicted by the Lindhard model 52 . This in¬ 
dicates that the computation of the least screened eigenpotentials might be avoided and carried out 
using model functions. 

If instead of one wishes to diagonalize %. the potential V'" arising from the rearrangements of 
the charge density in response to the applied perturbation needs to be included in the definition of 
the perturbation V' 1 *" of Eq. (Eq. (44)> 61 - 74 - 75 . In a generalized KS scheme the V ^ 1 is given by: 

ItfJar) = (Afy + ( I - a)At;+AVi + CtAVkxx) (47) 


where a is the fraction of exact exchange (EXX ) that is admixed to the semi local exchange poten¬ 
tial. The linear variation of the Hartree potential is 

AV;,|y/„k<T>= I </rAn,(r ) | [ / ^.- vWr) (48) 

and those of the exchange and correlation potentials are given by the functional derivatives: 




dV, 


x/c 


dn 


A/i.MvWr). 


mu 


(49) 
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The linear variation of the exact exchange potential (EXX) is expressed in terms of variations of 
the single particle orbitals 


= - I * /r/ [A^; k - CT <r'lV'.„k'rr(r)-E V'; rk - CT (r , )Av/,^ <T (r)] _____ V '»k«T< r’). 


HZ. 


(50) 


We note that calculations including V'/ 11 require a double self-consistent procedure (see Fig. Fig¬ 
ure 2): hence it is computationally more efficient to iteratively diagonalize x" first and then obtain 
the reducible polarizabiity X by linear algebra operations. We recall that both % a nd ar ° 
Hermitian operators ! 1 and because of Eq. (Eq. (18)) they have the same eigenvectors. 


Iiutinl guen 




Figure 2: (Color online) Diagrams of the iterative diagonalization of the irreducible (£°. left panel) 
and reducible <X- right panel) polarizability adopted in this work. In both cases the initial set of 
vectors {|p,)} are assigned with random components. At each iteration the polarizability matrix is 
computed by evaluating the density response to the i-th perturbation using Eq.(Eq. (44l)-(Eq. (46)). 
The two diagrams differ only by the self-consistent inclusion of V/* r in the solution of the Stern- 
heimer equation. 















































Figure 3: (Color online) Isosurfaces of (he square modulus of (he most screened eigenvector of the 
polarizability matrix of the silane molecule ( 0,(r)|’ in Fig. Figure 2, left panel). The iterative di¬ 
agonal ization was started from vectors with random components (left panel) and converges rapidly 
(3-4 iterations) to the potential shown in the right panel). 
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Figure 4: (Color online) Scalability of the PDEP iterative diagonalization (see Fig. Figure 2) of the 
static dielectric matrix of the COOH—Si/H,0 solid/liquid interface discussed in Sec. Section 4.2. 
The unit cell includes 492 atoms and 1560 valence electrons. 
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Figure 5: (Color online) Eigenvalues (A,) of the polarizability of H,0. SiH,. C h H 6 and N'aCI 
molecules, as obtained using the iterative diagonalization described in the left panel of Fig. Figure 2 
(see text), and adopting live different exchange-correlation potentials for the KS Hamiltonian (see 
Table Table 2). .V,,,,- is the number of valence bands 5 ’. 
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2.4 The evaluation of G and IV without empty electronic states using a Lanc- 
zos algorithm 

The calculation of the correlation self-energy. Eq.s (Eq. (27))-(Eq. (30)). and of the screening. 
Eq.s (Eq. (37))-(Eq. (40)). requires the computation of the matrix elements of defined in 

Eq. (Eq. (25)). for multiple values of 0). For each frequency 6). given two generic vectors | L) and 
l*> . we define 

*&(«) = (51) 

and 

«&(«>> = (L\/*°d° s (<o)f*°\R). (52) 


Eq. (Eq. (51)) can be easily written in terms of the eigenstates \y„ko and eigenvalues of the 
KS Hamiltonian: 


N*. 


"£*(») = -1 


(L\Y»ko)(Vnk<J \R) 


1=1 


irko 


- (0 


(53) 


Eq. (Eq. ( 52)) may be cast as well in terms of the occupied states and energies, by using the relation 
- 1 - F^ a and writing H$ s - P^ G H^ s F^ a (called deflated Hamiltonian) 


- X \R) 


(54) 


The Lanczos alorithm 7 s is used to obtain a set of Ni aneots orthonormal vectors Q — {|</,) : i — 1. A} 
that are used to recast the deflated Hamiltonian in Eq. (Eq. (54)) into a tri-diagonal form: 


( a \ 02 
02 <*2 01 


T = Q'HksQ = 



■ 



V 0» On/ 


(55) 
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where 


«" = («7 n |%k). (56) 

fin —II i^KS ~ a ») kn> “ fin l<//i-l) II • (57) 

The calculation of the sequence of vectors «/„). called Lanczos chain, is started by imposing |</|) — 
R): iterations are performed 79 by enforcing orthogonality through the recursive relation: 

*7n+l) = 77- [<«S “ U n) Vh<) ~ fin If/n-l). • (58) 

i w— 1 

The tridiagonal matrix T can be diagonalized using an orthogonal transformation U . so that l) — 
U'TU. Using to indicate the //-th element of the diagonal matrix D. we obtain: 

N/a «.<v I 

*K*(®) = - E {««, |R) ■ (59) 

.,i = l. "»2 10 

"2=1- 
' 1.1 = 1 

Because of the orthogonality of the elements belonging to a Lanczos chain, we have {r/„, | R) = 
4,i. yielding 

i 

= - E 1601 

0’ = 1 

Eq. (Eq. (60)) is written in a form similar to Eq. (Eq. (53)). where the coefficients of the expansion 
are independent of the value of (o and therefore it is not necessary to recompute them for multiple 
frequencies. However, the coefficients of the expansion in Eq. (Eq. (60)) depend by construction 
on the vector R) that is used to start the Lanczos chain. Therefore to evaluate \<o) one needs 
to solve as many Lanczos problems as the number of vectors R), while the eigendecoposidon used 
for in Eq. (Eq. (53)) is unique. Because Lanczos chains are independent of each other, 

the iterations can be performed in an embarrassingly parallel fashion, similarly to the procedure 
discussed in Sec. Section 2.3 for the computation of the PDEP basis set. 

In our calculations we used Eq. (Eq. (53)). with an explicit summation over occupied eigenstates. 
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for the evaluation of the terms in Eq. (Eq. (27)) and Eq. (Eq. (29)), whereas we used the Lane- 
zos expansion of Eq. (Eq. (60)) to evaluate the terms in Eq.s (Eq. (28)), |Eq. (30)), (Eq. (37))- 
(Eq. (40)). 


2.5 The contour deformation technique 

In Eq.s (Eq. (26))-<Eq. (30)) the frequency integration may be carried out by using complex analy¬ 
sis. and thus avoiding the integration in the real frequency domain. A closed integration contour on 
the complex plane is identilied for which Cauchy’s theorem and Jordan's Lemma apply 80 . This ap¬ 
proach is called contour deformation technique 29 81 and establishes a formal identity between the 
quantities reported in Eq. (Eq. (26)) and an equivalent set of quantities that are numerically more 
stable to compute. The poles of the Green’s function G£ s (r,r';&) + (o') are located at complex 
frequencies Z% ka . satisfying the relation 

-«'k<T = ~ <*> ~ /qsign(e„ ka - £ F ) (61) 

with a numerical residue given by 

Res {c^r.rVS**} = V'.*<r(r)i/4 0 <r'). (62) 

The poles of W p correspond to the plasmon energies of the system: cJJ — ±{£2 f , — iij). The matrix 
elements of the correlation self-energy can be computed by using the integration contour shown in 
Fig. Figure 6. yielding: 

-Moo 

l£(r.r':G>) = i j ^-G^.(r.r';d> + n/jM^r.r': o') + (63) 

—joo 

- L %k 0 (r)Vnk a (r)W p (rJ;£ ka )+ (64) 

+ Y. W*® V) W p (r, r' ; z ° ka ) . (65) 
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Figure 6: (Color online) Contours used in this work (see text). The integration contours I" and 
r enclose only the poles of the Green's function (dots) and exclude the poles of the screened 
Coulomb interaction ’J) (crosses). 

In view of the chosen contour, as the frequency to is varied, the poles of W p never fall inside the 
two closed contours T~ and T~. which therefore may only enclose poles of the Green's function. 
The correlation self-energy ls thus obtained as the sum of: i) an integral along the imaginary axis, 
where both G% s and W#pa are smooth functioas. and ii) all the numerical residues arising from 
<5*. shifted inside the integration contours, depending on the value of to. The matrix element of 
the correlation self-energy becomes: 

(Vntol ^c( E nk a )Wnka) = ~j ^ (V»kalG% s (r. E?£ + ito')W p (r.r'.ito') |% ka ) +(66) 

~ I fSZ ^(rlVV'dr.r':^,^ - £<&)*;*,(V) 

m 
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where the function is 


if e F < £ mka < E$f a 

if EZ < e mka < e F (67) 

otherwise 

Eq. (Eq. (67)) shows that only a finite number of residues needs to be computed 82 . Inserting 
Eq. (Eq. (22)) for \V p into Eq. (Figure 6). Eq.s (Eq. (26))-(Eq. (30)) are solved. The integration 
over the immaginary axis is performed numerically by considering a non-uniform grid, finer for 
small frequencies. With the introduction of the contour deformation technique we avoid the use 
of plasmon pole models 281 '’ 8 ' 86 to describe the frequency dependence of the screening and we 
overcome the limitations of the analytic continuation 8 ' 88 reported in Ref. [39.401. 

3 Verification and validation of results 

In this section we present several results obtained with the Go Wo method presented in Sec. Sec¬ 
tion 2. In particular we computed the vertical ionization potential (VIP) of closed and open-shell 
molecules and the electronic structure of several crystalline, amorphous and liquid systems. All re¬ 
sults were obtained by computing KS eigenvalues and eigenvectors with the Quant umEspresso 
package 8 *' and the Go Wo quasiparticle energies with the West code, which features a parallel im¬ 
plementation of the method of Sec. Section 2. 

3.1 Vertical ionization potentials of molecules 

We considered a subset of the G2/97 test set 90 composed of 36 closed shell molecules, listed in Ta¬ 
ble Table 3. Subsets of the G2/97 set were recently used to benchmark Go Wo calculations with both 
localized" 1 and plane wave 10 " 1 ' basis sets. Molecular geometries were taken from the NIST 
computational chemistry database" 7 , and no additional structural relaxations were performed. In 
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our calculations, each molecule was placed in a periodically repeated cell of edge 30bohr; the 
interaction between ionic cores and valence electrons was described by a PBE norm conserving 
pseudopotential: we used a plane wave basis set with a kinetic energy cutoff of 85 Ry (chosen so 
as to be appropriate for the hardest pseudopotential, i.e. those of oxygen and fluorine). At the 
DFT-KS level of theory we approximated the VIP by the absolute value of the highest occupied 
KS eigenvalue (HOMO)"' and we considered live different exchange and correlation functionals: 
PBE. PBEO. EXXc. B3LYP and HSE. whose expressions are summarized in Table Table 2. The 
computed DFT-KS VIP are reported in Table Table 3 within parentheses and in Fig. Figure 8 as 
crosses. As expected, hybrid functionals yielded a better agreement with experiments than PBE: 
the mean absolute relative errors (MARE) are 13.00%. 24.70%. 25.51% and 29.22% for EXXc. 
B3LYP. PBEO and HSE respectively, whereas the MARE of PBE is substantially higher. 37.98%. 
Corrections to the DFT eigenvalues were computed within the Go Wo approximation using the 5 
different starting points obtained with the various functionals. The PDEP basis set of each system 
was generated including a number of eigenpotentials Npj ep proportional to the number of valence 
electrons, for instance N p dt P = 1050 for the largest molecule considered here. i.e. C 6 H 6 . The 
(7(>Wo corrected VIPs are reported in Table Table 3 and in Fig. Figure 8 as dots: we obtained values 
in much better agreement with experiments, with MARE of 1.78%. 1.96%. 2.03%. 3.96% and 
4.49% for PBEO. B3LYP. HSE. PBE and EXXc starting points, respectively. We note that the QP 
corrections to HOMO DFT eigenvalues have different signs, depending on the starting point: the 
corrections lead to a decrease of the VIPs obtained with EXXc but to an increase of those computed 
with the other functionals. In Fig. Figure 7 we analyzed separately the matrix elements of V' lf . £y 
and E(- (see Eq. (Eq. (3))): the latter is the most affected by the choice of the exchange correlation 
functional at the DFT level. The matrix elements of L\ (panel <i) are instead weakly affected by 
the choice of the starting point. 

In many papers appeared in the literature 9, - 9S . Eq. (Eq. (3)) is solved using a linear approxima¬ 
tion 47 : 

Enka ~ e »'k<r +2«ko [(V'nko! ^(^nkir) I %ka) ~ {Ynko V'f.’ko)] (68) 
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where Z'. 1 — 1 — {\l/„ka £ a {(0) \ %ko) . Here we employed instead a secant method to 

find the roots of Eq. (Eq. (3)). where Eq. (Eq. (68)) was used to determine the starting point of the 
iterative procedure. The difference between VIPs obtained with Eq. (Eq. (68)) and using the secant 
method varies within 0 — 0.5 eV. see Fig. Figure 9. 

We also considered live open shell molecules, including the O, molecule in its triplet ground state. 
The VIPs computed at the DFT level using LDA or the PBE exchange-correlation functional", 
and by computing the QP corrections with GoWo are summarized in Table Table 4. The Go Wo 
results are in satisfactory agreement with the experimental data 97 and. for the systems considered 
here. LDA seems to provide a better starting point for Go Wo than PBE. 

We conclude this section dedicated to the validation of the West code for molecular sys¬ 
tems by showing that GoVV'o corrections may also improve higher order VIPs. i.e. vertical ion¬ 
ization energies obtained by removing electrons from single particle states deeper in energy than 
the HOMO. As an example we chose the thiophene (C,H,S) molecule whose spectral function 
A (to) — |^Tr {ImGo(to)} | was computed within GoWo. starting from DFT energies obtained using 
the PBEO functional (see Fig. Figure 10). We found that GoWo gives a much improved description 
of higher order VIPs with respect to KS-DFT. While for the experimental and the PBEO spectral 
functions we used an artificial smearing parameter of rj = 0.01 eV to simulate finite lifetimes, in 
the case of the GoWo spectral function the electronic lifetimes were computed from first principles, 
as the imaginary part of the electron self-energy. Our results are in good agreement with those 
reported by F. Caruso el al . 100 using localized basis sets. 
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Table 2: Exchange and correlation potentials used in this work (see text). For HSE. the screening 
parameter u — 0.106bohr 1 divides the exchange (x) contributions into short range (SR) and long 
range (LR) 101 . 


functional 

semilocal exchange 

nonlocal exchange 

corirlation 

PBE” 

V, PBC (r) 

- 

KT“( r) 

PBEO* 

0.75 <r) 

„.25Vf v *(r.r') 

V/“(r) 

EXXc* 

- 

vflr.r 1 ) 

v r«,r) 

B3LYP' 

0.08 V^“(r) + 0.72 V, m (r) 

0 . 2 V , l t '‘ v (r,r') 

0.19V.' XM (r)+0.8IVf"(r) 

HSE*' 

0.75V,' , ““(r:M)+V, r “^(r./i) 

„.25t' m '"(r.r’;p) 

V«*(r) 


" Ret. |102|. " Ref. [1()3|. 1 Kef. [KM]. 1 ' Kef. [105].' 

Ref. |10I | 

0 

■ PBE « PBE0 

■ EXXc ■ B3LYP 

■ HSE 




OP 

Figure 7: (Color online) The matrix elements of I',,.. L, and L,(£,, ko ) evaluated on the 
HOMO eigenstate, for different choices of the exchange and correlation potential (see Tab. Ta¬ 
ble 2). The bottom panel reports the QP correction, i.e. the difference (see 

Eq.s (Eq. (3)). (Eq. (7)>-<Eq. 110))). 
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Table 3: Vertical ionization potential (VIP, eV) of closed shell molecules. Experimental values are 
taken from the NIST computational chemistry database 97 . Each column reports the VIP obtained 
with the West code by performing CoWi calculations starting from the solutions of the Kohn- 
Sham equations with the exchange and correlation potential (see Tab. Table 2). specified within 
parentheses on the first row. In parentheses we report the absolute value of the HOMO energy 
prior to the application of (/oM'o corrections. ME. MAE. MRE and MARE stand for mean, mean 
absolute, mean relative and mean relative absolute error as compared to the experiment, respec¬ 
tively. 


Molecule 

r;„U’ (l (PBEl 

0„W'o(PBEO) 

O'dlVolEXXc) 

G„»H,(B3LYP) 

<; 0 U' 0 (HSE) 

Exp. 

C’H, 

11.10(7.20) 

11.38(8.43) 

11.66(12.19) 

11.37(8.45) 

11.30(8.03) 

11.49 

c 2 h 4 

10.35 (6.74) 

10.56(7.86) 

10.74(11.26) 

10.58(7.88) 

10.50(7.46) 

10.68 

c 4 h 4 s 

8.90(5.98) 

9.15(7.01) 

9.44 00.12) 

9.16(7.05) 

9.08 (6.62) 

8.86 

c„h 6 

9.10(6.33) 

9.32 (7.30) 

9.54 00.21) 

9.34 (7.33) 

9.25(6.91) 

9.25 

CH.Cl 

11.27 (7.10) 

11.57(8.50) 

11.89(12.81) 

11.56(8.57) 

11.49(8.11) 

11.29 

CH.OH 

10.47 (6.24) 

10.93(7.91) 

11.52 03.05) 

10.86(8.01) 

10.82(7.51) 

10.96 

CHjSH 

9.31 (5.55) 

9.57 (6.78) 

9.83 00.58) 

9.59 (6.87) 

9.49 (6.39) 

9.44 

ch 4 

13.99 (9.46) 

14.34 00.99) 

14.78 05.71) 

14.34 01.08) 

14.26 00.60) 

14.40 

ci 2 

11.48 (7.28) 

11.78(8.69) 

12.14(13.02) 

11.77(8.77) 

11.70(8.29) 

11.49 

CIF 

12.47 (7.83) 

12.84(9.441 

13.35 (14.33) 

12.83(9.53) 

12.76(8.04) 

12.77 

CO 

13.45 (9.06) 

14.01 (10.74) 

14.88 0 5.91) 

13.99 (10.88) 

13.91 (10.34) 

14.01 

CO, 

13.31 (9.08) 

13.73(10.69) 

14.34(15.77) 

13.65 00.76) 

13.64 00.29) 

13.78 

cs 

10.92 (7.38) 

11.53(8.89) 

12.51 (13.55) 

11.49(9.02) 

11.40(8.49) 

11.33" 

F, 

14.90 (9.42) 

15.51 (11.73) 

16.34(18.87) 

15.42 01.82) 

15.40 01.33) 

15.70 

HjCO 

10.38 (6.25) 

10.85(7.84) 

11.43 02.82) 

10.78(7.97) 

10.74(7.44) 

10.89 

HjO 

11.81 (7.23) 

12.37(9.04) 

12.91 (14.67) 

12.31 (9.12) 

12.24(8.63) 

12.62" 

HjO, 

10.96 (6.43) 

11.47(8.29) 

12.13 (14.06) 

11.40(8.40) 

11.36(7.88) 

11.70 

HCI 

12.54 (8.03) 

12.84(9.48) 

13.12 03.93) 

12.85(9.54) 

12.76(9.08) 

12.74- 

HCN 

13.30(9.02) 

13.63 00.39) 

13.96 04.54) 

13.60 00.40) 

13.55 (9.99) 

13.71 

HF 

15.14 (9.64) 

15.72 01.80) 

16.28 0 8.52) 

15.65 01.86) 

15.60 01.39) 

16.12 

HOCI 

10.93 (6.61) 

11.32(8.18) 

11.84(12.97) 

11.28 <8.291 

11.23(7.79) 

11.12- 

Li, 

5.03(3.231 

5.29 (3.80) 

5.39 (5.55) 

5.29 (3.87) 

5.14(3.43) 

5.11“ 

LiF 

9.97 (6.07) 

10.85(7.88) 

11.45 (13.74) 

10.79(7.95) 

10.66(7.47) 

11.30' 

LiH 

6.58 (4.35) 

7.57 (5.42) 

8.29 (8.86) 

7.51 (5.53) 

7.26 (5.02) 

7.90“ 

N, 

14.95 0 0.29) 

15.54 02.20) 

17.23 0 7.80) 

15.48 02.34) 

15.43 01.80) 

15.58 

n,h 4 

9.28 (5.28) 

9.72 (6.80) 

10.24(11.55) 

9.68 (6.92) 

9.62 (6.40) 


Na, 

4.73(3.13) 

4.86 (3.60) 

4.88 (5.04) 

4.89(3.72) 

4.78 (3.24) 

4.89" 

NaCI 

8.30(5.23) 

9.12(6.48) 

9.49 (10.47) 

9.09 (6.53) 

8.93(6.08) 

9.80 

NHj 

10.20(6.16) 

10.72(7.71) 

11.26 02.55) 

10.68(7.80) 

10.59(7.31) 

10.82 

P, 

10.44 (7.11) 

10.62(8.09) 

10.76 01.06) 

10.63(8.12) 

10.56(7.70) 

10.62 

PH ( 

10.46 (6.73) 

10.70(7.88) 

10.94(11.45) 

10.73(7.99) 

10.63(7.49) 

10.59 

SH, 

10.26 (6.29) 

10.53(7.55) 

10.76 01.40) 

10.55 (7.62) 

10.45(7.15) 

10.50 

Si,H 6 

10.45 (7.18) 

10.71 <8.281 

11.06 01.75) 

10.78(8.401 

10.64 (7.90) 

10.53 

SiH 4 

12.44 (8.52) 

12.82(9.86) 

13.32 0 4.03) 

12.83(9.97) 

12.72(9.46) 

12.30 

SiO 

11.09(7.49) 

11.51 (8.84) 

12.10(12.75) 

11.47(8.94) 

11.41 (8.441 

11.49 

SO, 

11.96(8.08) 

12.44(9.61) 

13.15 (14.39) 

12.39(9.75) 

12.34(9.22) 

12.50 

ME(cV) 

-0.42 (-4.29) 


0.49(1.50) 


-0.101-3.27) 


MAE (cV> 

0.44 (4.29) 

0.19(2.87) 

0.51 (1.50) 

0.21 (2.78) 

0.22 (3.27) 

- 

MRE(%) 

-3.68 (-37.98) 

0.15 (-25.51) 

4.31 (13.00) 

-0.02 (-24.70) 

-0.86 (-29.22) 

- 

MARE«‘J| 

3.96(37.98) 

1.78(25.51) 

4.49 0 3.00) 

1.96(24.70) 

2.03 (29.22) 

- 


" The NIST computational chcmisiry database* 1, docs not report the VIP but the ioni/ation potential. 


29 














Figure 8: (Color online) Comparison between calculated and experimental vertical ionizution po¬ 
tential (VIP) for the set of 36 closed-shell molecules listed in Tab. Table 3. Dots (crosses) refer to 
VIPs obtained at die Co Wo (DPT) level of theory. 



Figure 9: (Color online) Difference between the solution of Eq. (Eq. (3)) using a secant algorithm 
and employing the first order Taylor expansion of Eq. (Eq. (68)). 
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Table 4: Vertical ionization potential (VIP. eV) of open shell molecules. Experimental values are 
taken from the NIST computational chemistry database'' 7 . Each column reports the VIP obtained 
with the West code by performing Go l Vo calculations starting from the solutions of the Kohn-Sham 
equations with the exchange and correlation potential (LDA or PBE). specified w ithin parentheses 
on the first row. In parentheses we report the absolute value of the HOMO energy prior to the 
application of GqW'q corrections. 


Molecule 

spin 

GtWriLDA) 

<; ll n'„tPBE) 

Evp. 

CF 

0.5 

8.92(4.68) 

8.69(4.72) 

9.55 

NF 

1.0 

12.18(7.14) 

11.81 (7.051 

12.63 

NO? 

0.5 

10.82 (6.63) 

10.46(6.551 

11.23 

O? 

1.0 

12.11 (6.92) 

11.67(6.87) 

12.33 

s? 

1.0 

9.53(5.86) 

9.34 ( 5.82) 

9.55 
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Figure 10: (Color online) The spectral function A((ti) (see text) for the thiophene molecule 
(C. s H ,S). The peaks reported in the middle panel are located at the measured ionization poten¬ 
tials 1116 . The top (bottom) panel shows the spectral functions obtained at the <7 qM'o (DFT) level of 
theory, using the PBE0 exchange and correlation potential. The width of the peaks is set equal to 
0.01 cV. except for the top panel where electronic lifetimes are computed as imaginary part of the 
self-energies. 
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3.2 Electronic structure of crystalline, amorphous and aqueous systems 


We considered three crystalline systems Si. SiC and AlAs. one amorphous Si,N 4 . and one liq¬ 
uid water snapshot. The KS electronic structure was computed using super cells and the T point: 
64 atoms and 256 valence electrons for Si. SiC and AlAs, with cell edges of 20.53, 16.48 and 
21.34bohr, respectively; the amorphous Si ( N, sample consisted of 56 atoms and 256 electrons. 
For Si. SiC. AlAs and amorphous Si<N 4 we used a kinetic energy cutoff of 60Ry. The snapshot 
of 64 water molecules is taken from a 20 ps trajectory of a Born-Oppenheimer ub initio molecular 
dynamics simulation of liquid water (see Ref. ^41 ]): and it was described with a cutoff of 85 Ry. In 
our GoU'c, calculations for condensed systems we used — 1024. 

The QP energies of the crystalline solids at high symmetry k-points are reported in Table Ta¬ 
ble 5. Table 6 and Table 7. where KS energies are given within brackets. The results obtained with 
West compare well with those of other plane wave pseudopotential calculations and with experi¬ 
ments. 

The QP corrections of amorphous Si 3 N 4 and liquid water are reported in Fig. Figure 11. where 
again we found that the West results compare well with those of existing calculations' 11-10, and 
experiments 108 . 


Table 5: Quasiparticle (QP) energies of Si at high symmetry points, compared with previous cal¬ 
culations and experiment. 


k-point 

GoWolLDA) 

GyVoIPBE) 

rheo. 

Exp. 

L u 

2.26(1.47) 

2.29(1.59) 

2.2r. 2.14‘. 2.18*'. 2.13'. 2.19’. 2.05' 

2.1‘. 2.4±0.1" 

4, 

-1.251-1.20) 

-1.21 (-1.20) 

-1.23".-1.17*,-1.20*.-1.25'.-1.16* 

-1.2±0.2* 

r.5r 

3.35 (2.54) 

3.32 (2.55) 

3.25“. 3.24*. 3.24'. 3.23 d . 3.25'. 3.36'. 3.09* 

3.40*. 3.05' 


0.0 (0.0) 

0.0 (0.0) 

0.0 

0.0 

Xu 

1.44 (0.63) 

1.37 (0.72) 

1.36“. 1.41*. 1.34'. I.35 J , 1.31'. 1.43'. 1.01* 

1.3*. 1.25' 

Xi\ 

-2.92 (-2.87) 

-2.96 (-2.85) 

-2.88“. -2.80°. -2.80'. -2.83 d . -2.96'. -2.93'. -2.90 1 

-2.90 1 . -3.3±0.2“ 


" Ret. (40 ). b Ref. |551.‘ Ref. |87].‘ f Ref. |109|.' Ret. (104).' Ret. |l !()|.* Ret. |lll|. 
* Ref. (112). 1 Ret. (113|.' Ret. |114|. * Ref. |l 15 |. 1 Ref. [116|." Ref. |117) 
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Table 6: Quasiparticle <QP) energies of SiC ai high symmetry points, compared with previous 
calculations and experiment. 


k-poinl 

GuWi(LDA) 

GoWxPBE. 


Ev P . 

Li, 

6.46 <5.151 

6.37(5.19) 

6.43“. 6.30' . 6.45' 

6.35 7 

Li, 

-1.18 (-1.09) 

-1.161-1.06) 

-1.10“. -1.21* 

-1.15** 

Vu 

7.42 <6.29) 

7.52 (6.291 

7.26“. 7.19*. 7.23' 

7.4' 


0.0 (0.0) 

0.0 (0.0) 

0.0 

0.0 

X,. 

2.45 <1.29) 

2.28 (1-35) 

2.31“. 2.19*, 1.80' 

239*. 2.42* 

Xi, 

-3.46 (-3.25) 

-3.46 (-3.19) 

-3.47“. -3.53* 

-3.6* 


" Ref. H*>1- * Ref- |109).' Ref. |1 111.- 1 Ref. 11I2|,‘ Rrt. |1 ISj 


Table 7: Quasiparticle <QP) energies of AlAs ut high symmetry points, compared with previous 
calculations and experiment. 


k-poinl 

LDAi 

GbHbtPBE) 

Thro. 

Ecp. 

L\, 

3.08 (2.151 

2.94(2.15) 

3.02". 2.S4' 1 . 2.99 s 

2.36* 

Ly, 

-0.86 (-0.8(1) 

-0.90 (-0.84) 

-0.9". -0.87* 

- 

Vu 

3.15 (2.20) 

2.99 (2.23) 

2.96". 2.74°. 2.72 1 

3.13* 

Hlr 

0.0 (0.0) 

0.0 (0.0) 

0.0 

0.0 

Xu 

2.20(135) 

2.01 (1.37) 

2.13", 2.16*. 1.57 

2.23* 

x^ 

-2.25 (-2.15) 

-2.35 (-2.2 ll 

-230*. -2.27* 

-2.41* 


Ref. HOI. " Ref. |109|, * Ref. |1111. * Ref. (1121 “ Ref |I 19| 



Figure 11: (Color online) Quasiparticle (QP) corrections for a configuration of amorphous SijN 
(left panel) and liquid water at ambient conditions (right panel). 
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4 Large scale calculations 


The method discussed in Sec. Section 2. implemented in the West code and validated in Sec. Sec¬ 
tion 3 may be used to perform highly parallel Go Wo calculations and tackle large systems, with 
> 1000 of valence electrons in the unit cell. We discuss the performance of the method for both 
finite and periodic systems, in particular for Si nanocrystals and interfaces of functionalized Si 
surfaces and water, with up to 1344 and 1560 valence electrons in the unit cell, respectively. 

4.1 Silicon nanocrystals 

We considered four Si-NCs: Si 35 H, b (1.3nm>. Si 87 H 76 < 1.6nm). Si| 47 H |00 (1.9nm), and Si, 9 ,H )72 
(2.4nm) 120 . The structure of each NCs was obtained by carving out of bulk Si a sphere of Si atoms 
of given radius, by terminating all dangling bonds with H atoms and by relaxing the NC structure 
within DFT-PBE. A kinetic energy cutoff of 25 Ry. PBE norm-conserving pseudopotentials and a 
cubic cell of edge 90bohr were used. The computed HOMO and LU.VIO energies and the energy 
gap (Egjp) are reported in Table Table 8. The HOMO and LUMO energies referred to vacuum 
were obtained using the Makov-Payne 121 method. For each Si-NCs we considered — 2048. 
PDEP eigenvalues are reported in Fig. Figure 12 and they clearly show that the only difference 
between Si-NCs of different size appears for the most screened eigenpotentials. As discussed in 
Sec. Section 2.3, the PDEP eigenvalues of the least screened eigenpotentials are weakly affected 
by the microscopic structure of the system and may likely be predicted by model screening func¬ 
tions. The computed GqWq energy gaps for Si 15 H» (l . Si 87 H 7b . Si 147 H nf0 and Si, 93 H, 72 are 6.29. 
4.77. 4.20 and 3.46eV. respectively. These results are in good agreement with those of other re¬ 
cent calculations using MBPT or ASCF method. I22 . although our computed HOMO and LUMO 
energies differ slightly from those reported in Ref. [122]. 



Table 8: Quasipanicle (QP) energies and energy gap (E gap ) of Si nanocrystals. The Kohn-Sham 
eigenvalues obtained using the PBE exchange-correlation functional are reported in parentheses. 


Si-NC 

N« ( 

HOMO (eV) 

LUMO (eV) 

E» »V) 

GoW 0 (PBE) 

CqW 0 (PBE) 

G„W'„ (PBE) 

SijsH* 

176 

-7.59 (-6.08) 

-1.30 (-2J7) 

6.29(3.51) 

Si n H 76 

424 

-6.69 (-5.54) 

-1.92 (-2.96) 

4.77(2.58) 

r 

s 

£ 

688 

-6.48 (-5.44) 

-2.27 (-3.15) 

4.21 (2.29) 

Si.,,; II 172 

1344 

-5.82 (-5.19) 

-2.36 (-3.41) 

3.46(1.78) 
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4.2 Solid/liquid interfaces 


We now lurn lo discuss QP energies of extended, large systems. We considered two solid/liquid 
interfaces: H-Si/H,0 and C00H-Si/H,0. that were recently studied by T.A. Pham cl al . 121 
to align band edges of functionalized Sit 111) surfaces with water reduction and oxidation poten¬ 
tials. The orthorombic unit cell (L, x L v x L : ) of each system was obtained by interfacing 108 
water molecules with 72 Si atoms and by terminating the solid surface exposed to water with 24 
H atoms or 24 COOH groups, resulting in a (21.97 x 25.37 x 63.19)bohr , supercell with 1176 
valence electrons and a (21.97 x 25.37 x 67.53) bohr' supercell with 1560 valence electrons for 
H-Si/H,0 and for COOH—Si/H,0. respectively. Both interface geometries were extracted from a 
~ 30 ps trajectory of a Car-Parrinello molecular dynamics simulation of the interface where all wa¬ 
ter molecules and atoms of the semiconductor surfaces were allowed to move (see Ref. (123]) 124 . 
Side view's of the unit cells are shown in Fig. Figure 13, top panels. The KS electronic structure 
of both systems was obtained at the PBE level of theory using 85 Ry for the kinetic energy cutoff. 
The local density of states (LDOS) was obtained from the wavefunctions i/„ and energy levels £„ 
as 

LDOS( Zl £) = £|^|^|v'.( W )| 2 i(£-6,) (691 

where z is the axis perpendicular to the interface and 5 is the Dirac delta function. The LDOS 
of both systems, obtained at the PBE level of theory, is reported in Fig. Figure 13. middle panels. 
Those at Go Wo level, obtained by replacing the KS energies with QP energies in Eq. (Eq. (69)). are 
shown in Fig. Figure 13. bottom panels. The figures show that the method developed in Sec. Sec¬ 
tion 2 can be successfully used to compute the positions of the valence and conduction band edges 
of a realistic interface and hence to deline an electronic thickness of the interface, by analyzing 
how the bulk eigenvalues are modified in proximity of the interface. The method can of course 
be used for systems with impurity levels and to investigate semiconductor surfaces interfaced with 
aqueous solutions containing ions and to study the influence of ions on the electronic structure 
of the interface. The method developed here is not limited to solid/liquid interfaces and to pla- 
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nar geometries and has broad applicability to any complex < nanostructured > materials inclusive of 
heterogeneous interfaces 1 ". 



Figure 13: (Color online) The local density of states (LDOS. see text) of two solid/liquid interfaces: 
H-Si/H ,0 <left panels) and COOH -Si/H ,0 (right panels). The top panels report the side view of 
the unit cells. Bottom (middle) panels report the LDOS obtained using (7oW q (KS-DFT) energies 
in Eq. (Eq. (69)). A color scale that ranges from black to red is used to plot the LDOS: black areas 
indicate energy gap regions. 
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5 Conclusions 


We presented a formulation of the GW method for large scale calculations carried out with the 
plane-wave pseudopotential method. The evaluation of polarizabilities and electronic self-energies 
does not require the explicit computation of virtual states. Polarizabilities and dielectric matrices 
were represented with a basis set composed of the eigenstates of the dielectric matrix at zero fre¬ 
quency. obtained using iterative procedures. In the calculation of the correlation self-energy we 
avoided the use of the analytic continuation and carried out the frequency integration by means of a 
contour deformation technique. In addition we presented a parallel implementation of the method 
that allowed us to compute the electronic properties of large nanostructures and of solid/liquid in¬ 
terfaces. The method is not restricted to DFT inputs obtained with semi-local functionals but can 
be used in conjunction with DFT calculations with hybrid functionals. 

We presented a validation of the method for molecules (open and closed shell) and solids (both 
crystalline and amorphous) and found good agreement with data previously appeared in the liter¬ 
ature for converged calculations. We then applied our technique to silicon nanoparticles (up to a 
diameter of 2.4nm) and solid/liquid interfaces (with up to 1560 valence electrons in the unit cell). 
We showed that it is now possible to carry out many body perturbation theory calculation of real¬ 
istic slabs representing a semiconductor/water interface and to study in detail the modification of 
the bulk states at the interfaces and hence detine an electronic thickness of the interface. Work is in 
progress to couple our GW calculations with ab initio molecular dynamics simulations of realistic 
materials, and to include finite temperature and statistical effects in our MBPT calculations. 
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